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INTRODUCTION 


The  Gumbel  distribution  (extreme  value  type  I  for  maximumi  values)  [5;  Chap. 
3.51.  [6:  Chap.  21.4]  is  defined  by  its  probability  density  function  (PDF) 

f(x:X,b)  =  b"^  exp  -  (e"^  +  0  .  K  =  (x-X)b"'  ,  x  e  R  ,  (1.1) 

or  by  its  cumulative  distribution  function  (CDF) 

F(x:X.b)  =  exp  -  e"^  .  (1.2) 

with  parameters  shift  X  e  R,  and  scale  b  >  0.  It  is  of  considerable  importance  in 
various  statistical  data  interpretation  problems,  its  CDF  (1.2)  forms  the  basis  of  a 
simple,  highly  popular,  but  not  very  reliable,  linear  regression  procedure  for  the 
estimation  of  the  parameters  X  and  b  .  With  u  =  log  (-  log  F)  it  leads  to  the 
linear  function 

u  =  -b“'x  *  Xb“'. 

Therefore,  if  a  set  of  data  (xy,  Fy}  is  Gumbel  distributed,  the  points  (xy.  Uy  =  log 
(-log  Fy))  are  located  on  a  straight  line  in  the  (x.u)-plane.  In  other  words,  if  for  a 
given  data  set  {xy.  Fy),  the  points  (xy,  Uy)  seem  to  be  located  on  a  straight  line, 
the  unknown  parameters  of  the  Gumbel  PDF  (I.l)  can  be  determined  by  means  of  a 
least-squares  algorithm  [5:  Chap.  8.11. 

A  more  reliable  parameter  estimation  approach  can  be  based  on  the  maximum- 
likelihood  (ML)  principle  [14;  Chap.  12.51.  [8;  Chap.  5.4]  which  leads  to  a  simple 
algorithm  relative  to  the  Gumbel  PDF  (I.l).  It  is  suggested,  however,  to  apply  the  ML 
principle  to  a  more  general  PDF  class,  namely 

f(x;X,3,b)  =  &1.  oxp-Tfe^*^)  .  {;  =  (x-x)b''’  ,  xeR  .  (*') 

It  depends  on  the  parameters  shift  XeR.  shape  ^  >  0.  and  scale  b  >  0.  Introduction 
of  the  shape  parameter  3  makes  the  PDF  class  (**)  much  more  flexible  than  the 
original  Gumbel  ciass  (1.1)  which  is  contained  in  (**)  as  a  special  case  for  shape 
3  =  1. 


The  CDF  of  (**)  is  given  by 


r(x;X,iJ,b) 


f  (t;X,3,b)dt 


-oo 


where  Ffo.z)  is  the  incomplete  Gamma  function  with  lower  integration  limit  z  [3; 
6.350.21.  Cf  course,  this  CDF  leads  to  the  linear  relationship  explained  above  only  if 
3=1.  i.e.,  in  general,  a  simple  geometric  parameter  estimation  technique  is  no  longer 
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available.  It  may  be  of  interest  to  note,  however,  that  it  is  possible  to  show  that  log 
(1  -  r(bf;-X))  is  asymptotically  linear  in  t,  as  f  oo.  This  can  be  seen  by 
expressing  the  incomplete  Gamma  function  in  terms  of  the  degenerate  hypergeometric 
function  [3:  8.351.2]. 

It  is  the  objective  of  this  report  to  show  that  parameter  estimation  for  the 
three-parameter  generalized  Gumbel  class  (“)  via  the  ML  principle  can  be 
accomplished  by  means  of  a  reliable  algorithm  which  is  more  efficient  than  those 
currently  available  for  thr  classical  Gumbel  class  (1.1). 

Questions  of  hypothesis  justification  and  goodness  of  fit  are  outside  the 
framework  of  this  report  and.  thus,  are  not  addressed  here. 
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II.  DIFFUSION  CHARACTER  OF  THE  GENERALIZED  GUMBEL  DISTRIBUTION 
It  is  well  known  [91.  [10],  that  the  function 

g(x,t)  =  — j-2. - -  .  C  =  xc“’  ,  (II, ') 

X  >  0,  t  >  G,  U  >  0,  p  <  1 ,  with 


IS  the  Helta  function  initial  condition  (at  x  =  0.  t  =  0)  solution  of 
FokKor-Planck  equation 


(11.2) 


the  autonomous 


[a(x)  z(x,t)]^  -  [d(x)  z(x,t)]^  -  z^(x.t) 

A(x)  =  cxX^”'^  ,  (X  >  0 
D(x)  =  o<(2-0-p)x'  ^  -  rx 


0  ,  X  >  0,  t  >  0  . 

(diffusion  coefi  icient), 
(drift  coefficient)  . 


(11.3) 


With  t  =  (q  >  0  in  (11.2)  considered  as  a  parameter,  the  function  (11.1)  becomes  the 
lujporgamma  PDF  with  parameters  c  =  c  (to)  >  0.  p  <  1.  a  >  0  [101,  [1  ll. 

The  transformation  x  =  exp  -  [5;  (5.6)1  generates  from  (11.1)  the  function 

f(y,t)  =  g  (exp  -  2f'y,t)  I  dx/dy  I 

which,  after  y  has  been  replaced  by  x  .  is  of  the  form 

f  (x,t)  =  - ! -  exp-(c”^e  ^  ,  C  -  xb  '  ,  (lie) 

b  r  (^)  c^^ 


b  =  o  '2f  ,  3  =  (l-p)o  and  c  =  c(t)  given  in  (11.2).  For  any  fixed  t  -  to  >  0  and 
with  c  (to)  =  this  function  becomes  identical  with  (1.4)  (for  X  0). 


3 


The  transformation  x  =  exp  -  applied  to  the  differential  equation  (11.3) 
leads  (after  y  has  again  been  replaced  by  x  )  to  the  autonomous  Fokker-Planck 
equation 


[aV)  w(x,t)]^  -  [b^fx)  wCv.t)]^  -  w^(x.t)  =  0  . 
a“(x)  =  cxCJ^b^  exp  b"  ’  X  , 

D^fx)  =  cxa^b  (1-3)  exp  b”'x  +  axb  . 


of  which  the  function  (11.4)  is  a  solution.  By  means  of  different  parameter 
designations  the  function  (11.4)  can  also  be  transformed  so  that  it  becomes  a  nonlinear 
similarity  solution  of  (11.5)  [12:  Secs.  2,31. 

The  particular  case  1  -  p  =  a  reduces  the  PDF  (11.1)  to  the  Weibull  PDF  [n],  and 
(11.4)  becomes  the  Gumbel  PDF  with  3  =  1. 

111.  DATASETS 

It  will  be  assumed  that  the  given  statistical  data  are  of  the  general  form  of  a 
set  of  ordered  pairs  (x^.  fgv))  (v=l,  ...  M)  with  x^  <  X2  <  ...  <  Xfi  .  and  M  SN  = 

^^v=l  ^av  •  ^av  the  absolute  (integer  valued)  frequency  of  the  observation  xy. 
In  standard  statistical  terminology  data  given  in  this  particular  form  are  called 
histogram  data. 

It  will  furthermore  be  assumed  that,  for  histogram  estimation,  the  x^'s  are  the 
midpoints  of  class  intervals  [a  *  (v-l)Aa  .  a  *  vAa)  .  a  e  R,  Aa  >  0.  i.e.,  Xy  =  a  * 

(v-1 /2)Aa.  and  that  f0y  >0,  fgisl.  fgi-j^l.  Other  estimation  problems,  such  as 
popjiation  estimation  and  different  approaches  will  be  presented  elsewhere. 
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IV.  THE  LOGARITHMIC  LIKELIHOOD  FUNCTION 


For  the  generalized  Gumbel  PDF  C*)  with  parameter  vector  P  =  (X,3,d)  and  with 
the  abbreviation 

P;;(X)  =  Xy  ~  X 

the  likelihood  function  takes  the  particular  form 


(P)  -  3^b  ^(3)  exp-^ 


fg^exp-b  p, 


v-i 


The  f unct  ic  f-i 

<t(P)  -  n"'  logL(P) 


-  3  log  3  -  log  b  -  log  r(P)  -  3 


exp  -  b 


(IV.  1) 


where  fv  =  N''f3y  denotes  the  relative  frequencies  of  the  Xy's  with  E^y=ify=i,- 
is  the  logarithmic  likelihood  function  (LLF)  of  the  generalized  Gumbel  distribution 
class  (**). 

The  ML  principle  asserts  that  the  optimal  parameter  values  (if  they  exist) 
relative  tc  the  given  data  are  the  coordinates  of  the  point  P  =  (X,  S,  b)  in  the  cpen 
parameter  space  "P :  (X  e  R.  3  >  0.  b  >  0)  at  which  the  LLF  <|>(P)  takes  it  global 
maximum.. 
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V.  MAXIMUM-LIKELIHOOD  ESTIMATION 

In  this  report  ML  estimation  of  the  parameters  of  the  PDF  class  (**)  will  be  done 
by  means  of  the  derivative  equations  which  result  from  equating  to  zero  the  first 
partial  derivatives  of  the  LLF  <I>(P)  given  in  (IV.  1).  They  are  of  the  form 


f  = 


-i', 

V=1 


exp  -  b*'  p. 


=  0 


If.  =  log3-i 

op 


1  ^  ' 


1^  =  -b’’  *  3b'2  f^p^(i  -exp-b'’pv)  =  0  • 


where  'f(x)  =  d  log  r(x)/dx  is  the  psi  function  [3;  8.3621.  Their  appearance  can  be 
simplified  by  introduction  of  the  functions 


A(x.b) 


V=1 


exp  -  b“'  p  >  0  , 


B(A.b) 


V=1 


X  exp  -  b"’  p^  >  0 


V  V 


Then,  with  the  mean 


V=1 


^  x^  =  const,  .  0  <  X  <  Xj., 


they  take  the  form 


1  -  A  =  0  , 

log  p  -  1  -  -  b"'  (x  -  X)  -  A  =  0  . 

-b-3(x-X-B*XA)  ^  0  ,  (V.l) 
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By  means  of  the  first  one,  A  can  be  eliminated  from  the  other  two  so  that  only  the 
equations 


r(^)  -  b  ’  (x  -  X)  =  0 


-b  *  3(x-B)  =  0  .  (V.2) 

remain.  The  function 

r(3)  =  log  3  -  ^(3) 

which  appears  here  is  positive  for  3  ^  0  and  is  strictly  convex  and  decreasing,  r(3)  t 
00  as  3  i-  0,  r(3)  J-  0  as  3  t  oo,  as  can  be  verified  by  means  of  the  Laplace  integral 
representation  for  r(3)  [3;  8.361.81.  The  first  one  of  these  equations  allows  one  to 
express  b  in  terms  of  3  and  X, 


b-b(X.3)  =  (x-x)/r(3)>0  . 

Thus,  b  can  be  eliminated  from  equations  (V. I)  and  (V.2).  Setting 

A(x,b)  =  A(x,b(x,3)  =  Q(X,3)  =  £fyTy(X.3)  >  0  . 

v=l 

B(X,b)  =  B(X,b(X,3)  =  R(X,3)  =  £  f„T  (X.3)  >  0  . 


(V.3) 


v=l 


(V,4) 


(V.5) 


with 


T^  (X,3)  =  exp  - 


[ 


p^(x-X) 


(v=l . M) 


one  obtains 


g(X.3)  =  [x-R(X,3)]3r(3)  -  (x  -  X)  =  0 


(V.6) 


h(X,3)  =  1  -  Q(X,3)  =  0  . 


(V.7) 


If  (X,  3)  is  a  solution  of  this  system  of  equations,  then  b  follows  from  (V.3)  for  X 
=  X  ,  3  =  3.  The  objective,  therefore,  is  to  develop  an  efficient  numerical  algorithm 
for  the  solution  of  equations  (V.6)  and  (V.7).  Properties  of  the  functions  g(X,3)  and 
h(X,3)  which  are  essential  for  this  purpose  will  be  investigated  in  Section  VI. 
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Two  observations  relative  to  the  shift  parameter  X  are  of  importance.  Relation 
(V.3)  requires  that  X  <  x  for  the  scale  parameter  to  be  positive.  Furthermore,  if  X 

^  X)  =  min  {xp),  then  (x^  -  X)  (  x  -  X)"^  >  0  for  v  s  2  so  that  Q(X.3)  <  1  which 
contradicts  (V.7).  Consequently,  in  (V.6)  and  (V.7),  X  is  restricted  to  the  interval 

X|  <  X  <  x).  It  should  also  be  observed  that  the  positive  function  3r(3)  which 
appears  in  (V.6)  is  less  than  unity  for  3^0.  This  follows  from  the  integral 
representation  for  r(3)  mentioned  earlier. 


The  special  case  that  the  shape  parameter  3  is  assumed  to  be  unity  will  be 
discussed  next,  in  this  situation,  the  LLF  (IV.  1)  reduces  to 

<l>(X.b)  =  -logb-  ^f^(b~'p^  +  exp  -b"’p^]. 

v=l  '  * 


The  derivative  equations  now  lead  to 

1  -  A  =  0  . 

-b-B*  x=0  , 

with  A  and  B.  functions  of  x  and  b  .  defined  as  before.  The  first  equation  can  be 
solved  for  X  . 


X 


exp  -  b"'x^, 


v=i 

and  X  can  be  eliminated  from  B  so  that 


(V.8) 


B  = 


V=1 


■1-1 


^exp-b  x^ 


v=l 


f^x^exp-b  ’x.^, 


Therefore,  the  second  equation  becomes  an  equation  in  just  one  unknown. 


g(b)  =  X  -  b  - 


v=i 


exp-b''x^ 


v=i 


fyXyexp-b"'x^  =  0 


(V.9) 


It  is  easily  seen  that  g'(b)  is  negative  for  b  >  0  and  that  g(b)  Tx-X]  >0asbf0. 

and  g(b)  i  -  <»  as  b  t  <».  It  should  also  be  observed  that  (V.9)  implies  that  tS  <  x  if 
g(b)  =  0. 
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By  the  property  of  the  function  g(b)  just  discussed,  equation  (V.9)  has  exactly 
one  positive  root  b  .  It  immediately  gives  the  corresponding  value  X  for  the  shift 
parameter  from  (V.8)  for  b  =  1).  ML  estimation  for  the  two-parameter  Gumbel 
distribution  is.  therefore,  easily  accomplished  by  the  solution  of  the  single  equation 
(V.9). 
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VI.  PROPERTIES  OF  g(X,^)  AND  h(X.3) 

This  section  contains  a  description  of  properties  of  the  functions  g(X,3)  and 
h(X,3)  which  are  essential  for  the  design  of  an  efficient  solution  algorithm  for  the 
system  of  equations  (V.6)  and  (V.7).  The  function  h(X,3)  will  be  considered  first. 

By  means  of  the  substitution 

y  =  exp  -  (x-x)  r(3)  ,  (VI.i) 

which  takes  the  functions  Jy,  (X,3)  appearing  in  (V.4)  and  (V.5)  into 


T^(x.3)  =  u 


Xy-X 


equation  (V.7)  changes  into 


h(X,e) 


(VI.2) 


It  is  equivalent  to  the  equation 


x-x. 


t(y)  =  y 


y  =  u(y)  -  v(y)  =  0 


v=i 


For  fixed  X  in  the  interval  (x|,  x).  y  =  y(3)  becomes  a  one-to-one  mapping  of  (0,0°) 
into  (0.1).  Clearly,  f(y)  -*  0  as  y  t  i  (i.e.,  as  3  t  oo).  Furthermore, 


u(i)  ^  X  -  X,  >  0 


v’(i)  - 


£  rv(xv-xi)  =  X  -  X,  >  0. 
v=2 


Since  X  -  x]  >  X  -  X|  >  0.  it  follows  that  v' ( 1 )  >  u’ ( 1 )  >  0.  This  implies  that  0  < 
v(y)  <  u(y)  for  y  <  1,  y  sufficiently  close  to  1.  Consequently  f(y)  >  0  for  y  <  1.  y 
sufficiently  close  to  1,  and.  hence,  h(X,P)  >  0  for  3  sufficiently  large.  Furthermore. 
h(x,3)  0  as  3  °°- 
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Next.  u(y)  0,  v(y)  i  f  i  .  as  y  i  0.  Therefore.  f(y)  -  f  i  as  y  i  0.  and  hence. 

h(X,3)  i  -  00  as  3  i  0.  Consequently,  for  every  fixed  X  e  (X].  x).  the  function  h(X.3) 
has  an  odd  number  of  positive  zeros. 

Jt  can  now  be  shown  that,  for  every  fixed  X  e  (xj,  x).  equation  (V,7)  has  exactly 
one  positive  simple  root  ^  .  The  derivative  of  h(X.3)  with  respect  to  3  is  of  the 
form 

t  t 

(xQ-r)  =  -  q(3)  .  (V1.3) 

93  X  -  X  X  -  X 

R  and  Q  defined  in  (V.4)  and  (V.5).  respectively.  The  factor  -  (x  -  X)~^  r  is 
positive  for  every  3^0.  Suppose  the  second  factor 


q'3)  =  - 


£  fvT.,(xv-^) 
v=t 


has  two  distinct  positive  zeros.  0  <  3]  <  ^2-  By  Rolle’s  Theorem,  q'(3)  must  have  a 
zero  in  the  interval  (3],  32)-  But 


q’(3)  = 


0 


for  every  3  >  0.  Therefore,  q(3)  has  at  most  one  positive  zero.  By  properties  of  h 
discussed  earlier,  9h/93  has  at  least  one  positive  zero.  Therefore,  the  function 
q(3),  which  is  equivalent  to  9h/93.  has  exactly  one  positive  zero,  say.  3o-  Since  h 
has  an  odd  number  of  positive  zeros  and  since  h(X,3)  0  for  small  values  of  3  >  0 

and  h(X,3)  >  0  for  large  values  of  3  .  it  follows  that  h(X,3)  has  exactly  one 
positive  zero,  say,  3h(X),  such  that  0  <  3h(^)  3o-  2®'’°  ^h  's  simple  since 

q(3h)  0-  ii"  'iCPh)  =  0,  q(3)  would  have  two  distinct  zeros. 

It  is  now  possible  to  show  that  the  positive  function  3h(^).  unique  zero  of 

h(X,3)  for  given  X  e  (xj,  x  ),  implicitly  defined  by  equation  (V.7).  is  strictly 

monotonically  increasing  in  the  interval  (x^.  x  ).  Let  X  and  3  =  3h(^)  be  such 
that  h(X,  3)  =  0-  By  the  implicit  function  theorem  there  exists  a  function  3h(^) 
such  that  h(X,3h(X))  =  0  on  some  interval  (X).  X2)  which  contains  X  .  Since 
q(3h(^))  (^(3)  having  been  defined  in  (VI. 3))  is  different  from  zero,  in  fact 


q(3h(X))  =  XQ(X.  3h(^))  -  R(>^.  3h(>^))  =  3h(>^))  >  0  (V1.4) 
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for  every  X  e  (X|.  x  ),  the  function  3^^^)  <^3*^  continued  from  (X],  X2)  as  a 
differentiable  function  to  the  entire  interval  (x^.  x  ). 

Its  derivative  is  defined  by  the  identity 


Ah(x.wx))=il!!£^  [;-R(x.e„(M)] 
(x-x) 


x.3,(x))] 


dg^Cx) 

dx 


(V1.5) 


Since  S^v=l  (X.  3h(X))  -  1,  it  follows  from  Tchebychef’s  inequality  that  x  - 
R(X.  3h(X))  >  0  for  every  X  e  (xj,  x).  Furthermore,  by  (Vl.4),  X  -  R(X,  3h(X))  >  0. 
Therefore.  (VI. 5)  implies  that  d3h(X)/dX  >  0  for  X  t  (x],  x).  That  is  to  say,  3h(X)  is 
strictly  monotonically  increasing  in  the  interval  (x),  x). 

The  following  facts  concerning  the  function  3h(X)  should  finally  be  observed. 

By  (VI. 5),  its  derivative  is  given  by 


d3f,(X) 

r(VW) 

;;-r(x.  3f,(x)) 

■ _ 

1  ' 

1 

-x) 

The  positive  function  [  x  -  R(X,  3h(X))l  [X  -  R  (X,  3h(X))r'  J-  1  as  X  T  x  .  The 
positive  function  r(3h(X))  [-r  (3h(X))]~l  can  be  expressed  in  the  form 


r(x)  _  xr  (x) 

I  f 

-  r  (x)  x'f  (x)  -  1 


(Vi.7) 


Since  3h(X)  is  positive  and  strictly  monotonically  increasing,  the  argument  x  = 

3h(X)  in  (VI. 7)  cannot  go  to  zero  as  X  t  x  .  Furthermore,  for  x  >  0,  the  numerator 
xr(x)  in  (Vi.7)  is  bounded  above  by  unity,  and  xr(x)  J-  1/2  as  x  i  0°.  This  limit 
relation  can  be  established  by  means  of  the  Laplace  integral  representation  of  r(x) 
mentioned  before  and  by  application  of  Theorem  M.I.3  in  [2].  The  denominator 
function  in  tVl.7)  is  positive  for  x  >  0,  and  x'f'(x)  i  0  as  x  T  00.  This  limit  relation 
can  be  verified  by  use  of  an  integral  representation  of  'i'  (x)  [15:  Chap.  12.321. 
Consequently,  as  (VI. 6)  shows,  3V)(^)  T  00  as  X  t  x.  i.e.,  3h(^)  fias  a  singularity  at  X 

=  X  . 
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Now  suppose  would  approach  a  finite  limit  as  X  T  x.  Then,  as  (Vi.  1) 

shows,  y  i  0  as  X  t  x  .  i.e.,  the  second  term  in  (VI. 2)  will  not  remain  finite  as  X  T 
X.  This  establishes  a  contradiction  to  the  identity  (VI. 2)  with  3  =  3h(^)  3nd  X 
sufficiently  close  to  x  .  Therefore,  3fi(X)  T  oo  as  X  t  x  .  For  this  reason, 
numerical  solution  of  equation  (V.7)  may  fail  for  X  close  to  x  . 


Equation  (V.6)  will  now  be  considered.  The  transformation  (VI. 1)  changes  it  into 


g(x.3)  =  - 


-X  .  Xy 
X  -  y  2.  ^  x^  y 


P  (y)  log  y  -  1  =  0  . 


An  equivalent  equation  is 


-  X-Xi 
X  y 


P(y)  log  y 


x-x, 

U 


0 


(VI.8) 


in  which  X  is  considered  fixed  in  the  interval  (x,  x).  The  function  P(y)  cannot  be 
explicitly  expressed  in  terms  of  y  .  However,  it  is  only  necessary  to  know  that 


p(y)  log  y  =  -  ( x-x)  '  3r(3)  -*  { 


-(x-x)  as  3 -1  0,  i.e.,  as  y -1  0. 
1-1 


as  p  1 00.  j.e.,  as  y  T  1. 


The  second  limit  relation  has  essentially  been  established  already  in  connection  with 
(VI. 7).  The  first  one  follows  from  the  fact  that  xr(x)  T  i  as  x  i  0  [3:  8.362.2]. 

It  can  now  be  seen  that  the  left-hand  side  of  (VI.8)  approaches  positive  values 
as  y  t  0  and  as  y  T  i.  This  implies  that,  for  fixed  X  e  (xj,  x),  the  equation  g(X,3)  = 

0  has  an  even  number  of  positive  roots.  There  may  exist  values  of  X  e  (xj,  x)  for 
which  g(X,3)  =  0  has  no  positive  roots  at  all.  As  a  matter  of  fact,  a  more  careful 
investigation  of  equation  (VI.8)  reveals  that  there  may  exist  numbers  X|  and  X2  .  x; 

<  X]  <  X2  <  X  ,  such  that  g(X,3)  =  0  has  no  real  roots  for  X  e  (X],  X2). 
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VII.  THE  PARAMETER  ESTIMATION  ALGORITHM 

Estimation  of  the  parameters  of  the  generalized  Gumbel  distribution  class  C**), 

Section  I.  for  a  given  set  of  histogram  data  {(xy,  fy}  (v=l . M).  Xy  =  a  + 

(v  -  l/2)Aa,  is  based  on  the  solution  of  the  derivative  equations  (Section  V)  of  the 
associated  LLF  (IV.  1).  The  algorithm  presented  in  this  report  covers  the  three- 
parameter  unknown  case  as  well  as  the  two-parameter  unknown  case  in  which  the 
shape  parameter  3  is  assumed  to  be  equal  to  unity. 

In  the  three-parameter  case  the  equations  g(X.3)  =  0  (V.6),  and  h(X,3)  =  0 
(V.7)  must  be  solved.  Their  solution  (X.  $)  determines  the  scale  parameter  b  by 
(V.3).  In  the  two-parameter  case  the  single  equation  g(b)  =  0  (V.9)  must  be  solved. 
Its  solution  b  determines  the  shift  parameter  X  by  (V.8). 

1.  The  Three-Parameter  Case 

Suppose  that  (X.  is  a  solution  of  the  system  of  equations  (V.6)  and  (V.7) 

with  X  e  (x],  x)  and  ^  =  3h(^)  ^  0.  and  that  the  graphs  of  the  functions  implicitly 
defined  by  g  =  0  and  h  =  0  have  a  proper  intersection  at  (X,  3).  Then  there  are 

numbers  X|_  and  Xr  ,  X]  <  Xl  <  X  <  Lr  <  x  .  and  numbers  3l  =  3h(Xi_)  and  3r  = 
3h(XR)  .  0  <  3l  <  3r  .  such  that  h(XL.  3l)  =  0.  ^(Xr.  3r)  =  0  and  g(XL.  3l)  qUr.  3r) 

<  0.  Since  3h(X)  is  strictly  monotonically  increasing  (Section  VI),  the  solution 
point  (X.  3)  of  the  equations  (V.6)  and  (V.7)  is  boxed  in  the  interior  of  the  rectangle 
in  the  first  quadrant  of  the  (X.3)  plane  defined  by  the  points  (X|_,  3l).  (>^l>  3r)  • 

3l).  (^r.  3r). 

Boxing  of  the  solution  point  (X.  3)  requires  the  solution  of  the  equation  h(X,3)  = 
0  only  and  evaluation  of  the  sign  of  the  function  g(X,3). 

The  root  3l(^)  o*"  h(X,3)  =  0  is  first  bracketed  by  sign  change  detection  along 

the  search  sequence  3v  =  '0"’  (1.618)’^  (v=0,  1,  2 .  19).  Brent’s  method  [13: 

Chap.  7.3]  is  then  used  to  calculate  the  root.  In  the  numerical  algorithm  it  is  actually 
not  the  equation  h  =  0  which  is  solved  but  the  equivalent  equation  f  =  0  given  in  Sec. 
VI.  The  range  of  f  is  bounded  for  3  e  (O.oo)  whereas  the  range  of  h  is  not  bounded 
below  near  3  =  0.  The  use  of  f  instead  of  h  eliminates  possible  overflow  problems 
for  small  values  of  3  •  The  parameter  estimation  algorithm  is  aborted  if  no  root 
3h(X)  is  found  in  the  interval  (10“',  3 19)-  Shape  parameter  values  outside  this 
interval  would  lead  to  nearly  pathological  densities. 

To  obtain  X  values  for  boxing  the  solution  of  equations  (V.6)  and  (V.7),  the  3 
roots  of  h(X,3)  =  0  are  calculated  along  a  sequence  (X^®)  ,  X^0)r  ,  ,  x(')r  , 

X^'Y  .  -  )  of  X  values  with  X^^^  =  x  -  Aa/2  and 
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The  first  of  these  subsequences  converges  up  to  x  ,  the  second  down  to  xi  as  v  t  <». 
Let  _  3(v)p_  be  the  corresponding  roots  of  h  =  0. 

The  sign  of  g(X,3)  is  evaluated  at  the  points  3hd 

3^^Y).  The  process  is  terminated  if  at  two  successive  points,  either  in  the 
{x(^)r)  or  in  the  {x(''^)l}  subsequence,  the  product  of  the  g  functions  is  negative, 
i.e.,  if  boxing  of  (X,  3)  has  been  achieved. 

The  boxing  rectangle  is  subsequently  refined  by  bisection  of  the  X  interval 
until  the  accuracy  requirement  (AX  <  10"^)  n  (A3  <  10*^)  is  satisfied. 

Monotonicity  of  the  function  3h(^)  is  not  exploited  in  the  actual  numerical  algorithm 
because  of  the  interplay  between  the  one-dimensional  Brent’s  method  and  the  two- 
dimensional  boxing  process. 

2.  The  Two-Parameter  .Case 

The  root  b  of  equation  (V.9),  g(b)  =  0.  is  first  bracketed  by  evaluation  of 

g(bv)  along  the  search  sequence  by  =  10"'  (1.618)^  (v=0.1,2 .  19).  The  reason  for 

the  restriction  of  b  to  this  interval  is  analogous  to  that  given  earlier  for  3  . 

Remark.  For  computational  convenience  the  transformation  x^  -*  x^,  -  X]  is 
applied  before  the  start  of  the  actual  computations.  It  affects  only  the  shift 
parameter.  The  algorithm  calculates  the  relative  shift  value  Xq  and  returns  the  true 
shift  value  X  =  Xg  +  X] . 
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VIII. 


EXAMPLES 


To  demonstrate  the  algorithm  a  number  of  examples  are  presented. 

Accompanying  tables  are  given  in  Section  IX. 

DiamaJa,  1 

N  =  07  observations  of  annual  24-hour  maximum  rainfalls  (in  points)  at  Sidney, 
Australia,  over  the  period  1859-1945  are  given  in  Table  1.1  [7].  Grouping  into 
histogram  absolute  frequency  data  has  been  performed  on  the  data  of  Table  1.1  in  five 

different  ways,  Gy  ,  as  displayed  in  the  second  columns  of  Tables  1.3Gy  (v=I . 5). 

The  first  column  shows  the  class  interval  numbers  v  .  The  class  interval  data  are 
given  in  Table  1.2. 

Parameter  estimation  has  been  performed  with  three  and  with  two  (3  =  0 
parameters.  The  results  are  shown  in  Tables  1.4Gy  (v  =  1 . 5). 

The  estimated  parameter  values  are  used  to  calculate  the  expected  absolute 
frequencies  from  the  PDF  (**).  For  the  three-parameter  estimates  they  are  given  in 
the  third  columns  of  Tables  1.3Gy  ,  and  for  the  two-parameter  estimates  in  the  fifth 
columns.  Next  to  the  calculated  frequencies  in  Tables  1.3Gy  are  shown  the  chi- 
square  values.  As  mentioned  in  the  Introduction,  no  significance  will  be  attached  to 
them  within  the  framework  of  this  report. 

ExamBla,  2 

A  total  of  N  =  69  observations  of  frequencies  of  annual  maxima  of  rainfall  (in 
inch)  in  24  hours  at  Camden  Square,  London,  over  the  period  1860-1948  [1:  Chap.  8.7] 
have  been  grouped  into  M  =  12  class  intervals  (a  +  (v-l)Aa  .  a  +  vAa)  with  a  =  1/2, 

Aa  =  1/4.  The  data  are  given  in  Table  2.1.  The  results  of  the  three-parameter  and 
two-parameter  estimates  are  shown  in  Table  2.2.  The  calculated  expected 
frequencies  together  with  the  corresponding  chi-square  values  are  displayed  in  Table 
2.1. 


Example  3 

This  example  concerns  the  distribution  of  the  greatest  ages  of  men  [4:  Sec.  I, 
Table  I,  2nd  ed.l.  A  total  of  N  =  52  observations  have  been  grouped  into  M  =  9  class 
intervals  [a  +  (v-l)Aa,  a  +  vAa)  with  a  =  95.5  and  Aa  =  1.  Table  3.1  shows  the  given 
data.  Calculated  results  are  shown  in  Table  3.2  (parameters)  and  in  Table  3.1 
(frequencies  and  chi-square  values). 

Example  4 

The  last  example  deals  with  the  ensemble  distribution  of  the  greatest  ages  of 
men  and  women  [4:  Sec.  1,  Table  1,  4th  col.].  The  total  number  of  observations  is  N  = 
104  corresponding  to  M  =  II  class  intervals  (a  +  (v-l)Aa,  a  +  vAa),  a  =  95.5,  Aa  =  1. 
Table  4.1  shows  the  given  data.  The  results  are  presented  in  Tables  4.2  and  4.1. 
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IX.  TABLES 


17/(18  Blank) 


Table  1.1  Annual  haximum  Twenty-Four  Hour  Rainfalls  (in  Points)  at  Sidney. 
Australia.  1859-1945. 
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Tablel.2  Class  Intervals  and  Number  of  Classes  for  Groupings  Gy  (v=1 . 5) 
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Table  1.3G]  Grouped  Frequency  Data  for  Grouping  G] 


GUMBEL  DISTRIBUTION  CALCULATIONS 
Data  file  is  LNGDSET1 .DAT 
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2.46 

Total  s 

87 

86.43 

86 . 08 
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Table  1.3G2  Grouped  Frequency  Data  for  Grouping  G2 


GUMBEL  DISTRIBUTION  CALCULATIONS 
Data  file  Is  LN606ET2.0AT 


3  Para 

Estimate 

2  Para 

E  s  t ima  t  e 

V 

f (  abs) 

f(  cal) 

X**2 

f(  cal) 

X»»2 

1 

8 

6.79 

.22 

6.79 

.22 

2 

IB 

20 . 64 

.  34 

19.70 

.  15 

3 

25 

22.61 

.  25 

22 . 74 

.  23 

4 

16 

16.10 

.  00 

16.78 

.  04 

B 

7 

9.56 

.  68 

9.97 

.89 

6 

7 

5.25 

.58 

5.34 

.52 

7 

3 

2.79 

.  02 

2.71 

.03 

8 

2 

1 .47 

.  19 

1 . 35 

.32 

9 

0 

.  77 

.77 

.66 

.66 

10 

1 

.  40 

.91 

.  32 

1 .43 

Totals 

87 

86 . 38 

86 . 35 

Table  I.3G3  Grouped  Frequency  Data  for  Grouping  G3 


GUMBEl  DISTRIBUTION  CALCULATIONS 
Data  file  is  LNGDSET3.DAT 


1 

2 

3 

4 

5 

6 
7 

a 

9 

10 
1  1 
12 

13 

14 
16 
1  6 
1  7 
18 

19 

20 

Totals 


f(  aba) 
6 

5 

6 
15 
12 
1  1 
10 

0 

5 

5 

3 

1 

3 

1 

1 

0 

0 

0 

0 

1 

87 


1  Para 

Estimate 

2 

Para 

Estimate 

ca  J 

L) 

f ( cal) 

X*>2 

3, 

.  13 

2 . 63 

3, 

.20 

2.46 

7 

.  28 

.  7  1 

6 

,  74 

.45 

10 

.  80 

.  72 

10 

.  04 

.42 

12 

1  6 

.66 

1  1 

.  77 

.69 

1  1 

.  56 

.  02 

1  1 

.69 

.01 

9 

.  90 

.  12 

10 

.  36 

.  04 

7 

.  95 

.53 

8 

.49 

.  27 

e 

.  14 

6.14 

6 

.  59 

6 . 59 

4 

.  63 

.  03 

4 

.93 

.00 

3 

.  44 

.71 

3 

.59 

.55 

2 

.  53 

.09 

2 

.57 

.07 

1 

.  85 

.  39 

1 

.82 

.  37 

1 

.  35 

2.02 

1 

.28 

2 . 32 

.  98 

.  00 

.89 

.01 

.  7  1 

.  1  1 

.  62 

.  23 

.  52 

.  52 

.  43 

.43 

.  38 

.  38 

.  30 

.  30 

.  27 

.  27 

.21 

.21 

.  20 

.  20 

.  14 

.  14 

.  14 

5.11 

.  10 

8.34 

85 

.  93 

85 

.  74 
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Table  1.3G^}  Grouped  Frequency  Data  for  Grouping  G^ 


GUMBEL  DISTRIBUTION  CALCULATIONS 
Data  file  is  LNGDSET4.DAT 


3  Para 

Estimate 

2  Para 

Estimate 

V 

f (  abs) 

f( cal) 

X*«2 

f(  cal) 

1 

2 

1 .84 

.01 

2.08 

.  00 

2 

9 

8.91 

.00 

6. 55 

.02 

3 

IS 

16.12 

.08 

15.46 

.01 

4 

20 

17 . 39 

.39 

17.35 

.  40 

5 

13 

14 . 32 

.  12 

14.76 

.21 

6 

e 

10.23 

.40 

10.70 

.  68 

9 

7 

6.70 

.01 

7.06 

.  00 

s 

6 

4 . 33 

.  6B 

4.41 

.57 

9 

3 

2.70 

.03 

2 . 67 

.  04 

10 

2 

1 . 67 

.  06 

1 .59 

.  10 

1  1 

1 

1 . 03 

.00 

.94 

.  00 

1  2 

0 

.63 

.63 

.55 

.55 

1  3 

0 

.  39 

.39 

.  32 

.  32 

1  4 

1 

.  24 

2.47 

.  19 

3.56 

1  0 1  a  1  s 

87 

86 . 57 

86.62 

Table  I.3G5  Grouped  Frequency  Data  for  Grouping  G5 


GUMBEL  DISTRIBUTION  CALCULATIONS 


Data  file  is  LNGDSETB.DAT 

3  Pa  a  Es 


V 

f (  abs) 

f(  cal) 

1 

6 

3 . 98 

2 

8 

12.24 

3 

20 

17.51 

4 

17 

16.70 

5 

16 

12.02 

6 

4 

0.70 

7 

6 

5 . 66 

8 

4 

3.53 

9 

3 

2.16 

10 

1 

1.31 

1  1 

1 

.  79 

12 

0 

.48 

13 

0 

.29 

i4 

1 

.  17 

otal  8 

07 

86.43 

mate 

2  Para 

Estimate 

X«-M-2 

f(  cal) 

X*-»2 

1 . 02 

4.05 

.  93 

1 . 47 

11.76 

1 .20 

.  35 

17.16 

.47 

.01 

16.88 

.00 

.79 

13.23 

.50 

2.60 

9.09 

2.85 

.02 

5.80 

.01 

.06 

3.54 

.06 

.  32 

2.11 

.38 

.07 

1 .24 

.05 

.  05 

.  72 

.  1  1 

.46 

.42 

.42 

.29 

.24 

.24 

3.94 

.  14 

86 . 36 

5.33 

21 


Table  1.4Gi  Three-  and  Two-Parameter  Estimates  for  Grouping  G] 


Input  file  is  LNGOSET1.DAT 
This  is  a  3  parameter  fit 

#*■«■•«•« -ft**  *■»*•»■«■  •!»■■»*«•■«•■»»■»«* -tnnm- -it  -tn;- 

THE  FINAL  VALUES 

Lambda  >  .315S19E+ee3 

b  -  .762429E■^002 

Bata  -  .4369i6E-«-000 

innnj- -jnnnf  iHf- 

Input  file  is  LNGD6ET1.DAT 
This  Is  a  2  parameter  fit 

«■»*«**•»*«*»«■» -JHHf  ** -tf  ■» -tt  -tf  *■»■»* -tUf  -If  « -if  ■IHHf 

THE  FINAL  VALUES 

Lambda  -  .348179E+003 

b  -  .  134030E-«-003 

Beta  =  .  1 00000E-«-00 1 


Table  1.4G2  Three-  and  Two-Parameter  Estimates  for  Grouping  G2 


Input  file  is  LNGDSET2.DAT 
This  is  a  3  parameter  fit 

THE  FINAL  VALUES 

Lambda  =  .336618E+003 

b  -  .116041E+003 

Beta  =  .759879E+000 

-M- ■»■*■«•  tt -ii- •»*■«■**  -innm- ■»*■» -H- 

Input  file  is  LNGDSET2.DAT 
This  is  a  2  parameter  fit 

•it  * -inm- ii- -innnf  -iH}**  ij- -tnni- * -snnnnnm- * 

THE  FINAL  VALUES 

Lambda  -  .346930E+003 

b  -  .138075E+003 

Beta  =  .  100000E-^001 

if  it -it  -it « -if  if -it  •» -M- -ft  *■!}■■»•»•«■  ■it  tt -it  ** -it  it -M- * -it  •» -it 
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Table  I.4G3  Three-  and  Two-Parameter  Estimates  for  Grouping  G3 

Input  file  is  LNG0BET3.DAT 
This  is  a  3  parameter  fit 

THE  FINAL  VALUES 

Lambda  »  .331309E+003 

b  -  .102070E+003 

Beta  -  .661333E-t-000 

Input  file  ie  LNG0SET3.DAT 
This  is  a  2  parameter  fit 

THE  FINAL  VALUES 

Lambda  «=  .346724E  +  003 

b  =  .13413SE+003 

Beta  =  .  1 00000E-t-00 1 

il- *«■■»••»•»»•»•«■»•** -jj- -it  *»*■«■•«•**■«■  * -M-t;- 


Table  I.4G4  Three-  and  Two-Parameter  Estimates  for  Grouping  G4 


Input  file  ie  LNGD6ET4.DAT 
This  is  a  3  parameter  fit 

•»■»*•»  innni' •»■»•»*•»•*•»**■»•«■■»•»«■*«■■»•«•«•**■»  w -ft* -a- * 

THE  FINAL  VALUES 

Lambda  —  ■  339635E-*-003 

b  -  .116049E+003 

Beta  >  .762044E+000 

** -a- -a- -a  * 

Input  file  is  LNGDSET4.DAT 
This  is  a  2  parameter  fit 

a  * -(nt  * -jnnnnnnm- -a- *  « •a- ■»■  •»  » •»** -a- * -a- « -a  * -a 

THE  FINAL  VALUES 

Lambda  “  .349790E+003 

b  -  .137703E+003 

Beta  .  1 00000E-t-00 1 

■a  a  a- -a -a -a  a  a- -a -a  a- a- -a -a  a- a^  * -a -a  a  a -a -a -a  a- -a  a  a -a -a  a- a -a  <i- -a -a  a- a- -a  ■>; 
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Table  I.4G5  Three-  and  Two-Parameter  Estimates  for  Grouping  G5 

Input  file  is  LNQD6ETS.DAT 
This  is  a  3  parameter  fit 

THE  FINAL  VALUES 

Lambda  -  .336333Efe03 

b  -  .118707E+003 

Beta  -  .60SS87E«000 

************»•)»■•»***•»■»■***•»■»*•»«■«••»***■»****■#••«•■» 

Input  file  is  LNG0SET5.DAT 
This  Is  a  2  parameter  fit 

THE  FINAL  VALUES 

Lambda  ■  .3441S7E+003 

b  “  .13Sai1E+003 

Seta  «  .  10000eE-i-0ei 

*** -tt  **«■»■•»**  HJ- ■»■***  <nnnm' * -H- ■»**** -M- ■»**»* -tut 


Table  2.1  Annual  Maxima  of  Rainfall  (in  Inch)  in  24  Hours  at  Camden  Square,  London, 
1860-1948 

GUMBEL  OISTRIBUTIGN  CALCULATIONS 
Data  file  Is  GL0SET2.DAT 


3  Para 

Estimate 

2  Para 

Estimate 

V 

f  (  abs) 

f(  cal) 

X*»2 

f (  cal) 

X**2 

1 

4 

3.58 

.05 

3.86 

.00 

2 

18 

19.31 

.09 

18.34 

.01 

3 

24 

2S.08 

.05 

25.07 

.05 

4 

24 

18.43 

1 .68 

19.21 

1  .  19 

S 

10 

10.79 

.06 

11.26 

.  14 

6 

4 

6.77 

.54 

5.82 

.57 

7 

2 

2 . 97 

.32 

2.84 

.25 

6 

0 

1.51 

1.51 

1 . 35 

1 .35 

9 

0 

.  76 

.76 

.63 

.63 

10 

2 

.  38 

6.79 

.  30 

9.82 

1  1 

0 

.  19 

.  19 

.  14 

.  14 

12 

1 

.  10 

8.37 

.  06 

13.69 

ota  1  8 

89 

88 . 87 

88.89 

24 


Table  2.2  Three-  and  Two-Parameter  Estimates 


Input  file  is  gldaet2.dat 
This  is  a  3  parameter  fit 

THE  FINAL  VALUES 

Lambda  -  .  1  e7645E-t-00 1 

b  -  .272609E-^a0e 

Beta  °  .749054E-t-0e0 

ft*ftft***ftft***ft*ftft*ft*fttt*ftftftftftftftftft*ftftft«**ftftft 

Input  file  is  gldset2.dat 
This  is  a  2  parameter  fit 

ftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftft 

THE  FINAL  VALUES 

Lambda  >  .  110330E-t-001 

b  -  .  32563 1E-)-000 

Beta  -  .  100000E-t-e01 

ftftftftftftftftftftftftftftftftft-ftftftftftftftftftftftftftftftftftftftftftftftft 


Table  3. 1  Greatest  Ages  of  hen 

GUMBEL  DISTRIBUTION  CALCULATIONS 
Data  file  is  GLDSET1.0AT 


1 

2 

3 

4 

5 

6 
7 
6 
9 

Totals 


f  (  ab s) 
1 
e 

9 

9 

10 

5 

2 

7 

1 

52 


3  Para 

Estimate 

f ( cal) 

X*ft2 

1  .56 

.  20 

5.57 

1 .06 

9.00 

.  06 

10.93 

.  34 

9.  12 

.09 

6.31 

.27 

3.06 

.91 

2.21 

10.41 

1 .20 
50,57 

.03 

2  Para 

Estimate 

cal) 

X»»2 

1 .37 

.  10 

5.94 

.71 

10.38 

.  10 

10.93 

.34 

0.69 

.20 

6.90 

.  14 

3.67 

.76 

2. 10 

10.70 

1 .25 

.05 

50.31 
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Table  3.2  Three-  and  Two-Parameter  Estimates 


Input  file  is  6LD6ET1.0AT 
This  is  a  3  parameter  fit 

THE  FINAL  VALUES 

Lambda  ™  .  98808  IE-i-602 

b  -  .22424eE-i>e01 

Bata  -  .  1S6860E-I-OO1 

Input  file  is  GLD6ET1.DAT 
This  is  a  2  parameter  fit 

THE  FINAL  VALUES 

Lambda  -  .  9862 1 4£ •*-002 

b  -  .191052E+001 

Bata  «  .  10000eE-t-801 


Table  4.1  Greatest  Ensemble  Ages 


GUMBEL  DISTRIBUTION  CALCULATIONS 
Data  file  is  GLDSET3.DAT 


3  Para 

Estimate 

2  Para 

Estimate 

V 

f (  aba) 

f(  cal) 

X**2 

f(  cal) 

X«*2 

1 

1 

1 . 99 

.49 

1 .20 

.03 

2 

1  1 

7 . 39 

1 . 77 

6. 12 

1  .02 

3 

IS 

1S.S6 

.02 

18.31 

.60 

4 

18 

21.27 

.50 

22.29 

.83 

S 

22 

20.93 

.05 

19.23 

.40 

6 

IS 

16.04 

.07 

13.69 

.  12 

7 

10 

10.  16 

.00 

8.76 

.  18 

8 

6 

S.S7 

1 .06 

S.28 

1  .41 

9 

3 

2.73 

.03 

3.07 

.00 

10 

0 

1 .23 

1 .23 

1 .76 

1  .76 

1  1 

1 

.52 

.44 

.99 

.00 

Totals 

104 

103.39 

102.70 
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Table  4.2  Three-  and  Two-Parameter  Estimates 


Input  flla  it  6LD6ET3.DAT 
This  is  a  3  paramatar  fit 

THE  FIKAL  VALUES 

Lambda  -  .  393298E-i-0ei 

b  -  .366330E+eoi 

Bata  -  .  3896 1SE -^001 

Input  fila  is  GL0SET3.0AT 
This  is  a  2  paramatar  fit 

THE  FINAL  VALUES 

Lambda  -  .  347478E<<-001 

b  -  .171638E+e01 

Bata  -  .  100e0eE-i>001 


27/(28  Blank) 


X.  PROGRAM  LISTING 


29/(30  Blank) 


o  o  o 


C  Research  Directorate,  RDE.C,  USA  MICOM 

Or .  C  E  Hall,  Jr . 

1  Sept  1988 
mod  4  Jan  89 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-2J 

COMMON  /GGG/  NCL, FREL( 225) ,RHO( 225) ,XC( 225) ,SUMA, SUMS, XBAH 
COMMON  /SSS/  OSHIFT,  NTOT,  ALPHA,  J20R3 

CHARACTER*30  FILENAME 
CHARACT£H*?5  COMMENT 

C  Send  standard  output  to  the  printer 

0PEN(7,  FILE  =  'PRN') 

C  Open  i->put  file 

CALL  HISTOHDC FILENAME , UELTAA) 

C  Output  the  results  both  to  the  screen  and  to  the  output  tile 

WRITE(7,*)  '  Input  file  is  FILENAME 

C  I/O  for  a  two  or  three  parameter  fit 

WRITEf^-^,'fA\)')  '  Enter  desired  parameter  fit,  2  or  3-' 

READ(  ■•- ,  •■=•)  J2QR3 

WRITEC7,'(A,I2,A)')  '  This  is  a  ',J20H3,  '  parameter  fit' 

C  Calculate  XSAR 

CALL  CALX8AH 

IF  CJ20R3.EQ.3)  THEN 

i:  Calculate  the  roots  of  g(  .  )  and  hi.) 

lERR  =  a 

CALL  ROQTGC XRQQT , Y800T , lERR) 

IF  (  lERR . NE . 0)  GOTO  900 

C  Calculate  the  value  of  b 

TEMPI  »  OLOG( YROOT)  -  PSI(YROOT) 

TEMP2  =  XSAR  -  XROOT 
QVAL  =  TEMP2  /  TEMPI 
REAM  =  XROOT  -  OSHIFT 
ELSE 

C  Scale  BVAL  to  match  HOOTH  routine 

IF  (,  XCl  NCL)  .  GT  .  70 . 00)  THEN 

ALPHA  -  XC(,  NCL)  /  70.00 
ELSE 

ALPHA  =  1 . Du 

cnoiF 

t  Sol  -or  d '/ Al  ,  vi  j  HJUIH 

CALL  HOO  THi;  X  ,  BVAL  ,  lERR) 

IF  llERR.NE.0)  GOTO  200 

C  Calculate  RlAM  and  set  BETA  l  YROOT) 

BVAL  =  OVAL  ALPHA 

SUM1  =  0.00 

00  100  ICNT  =  1,  NCL 

bUMI  =  SUMl  +  FRELI  ICNT)  ,  0 1 XPl  XCC  I ON T )  /  8 V Al ) 

130  CONTINUE 


31 


c 

90a 

999 

C 


c 


1 00 

200 


RLAM  -  (-1.D0)  »  BVAL  *  OLOG(BUMI)  -  DSHIFT 
YROOT  =  1 . D0 
ENDIF 


Write  values  to  the  screen 


WRITEC  9 , 
WRITE! 7,») 
WRITEC  7, 
WRITEC  7,  'C 
WRITEC  7,  'C 
WHITE! 7,  'C 
WRITEC  7 , 


'  THE  FINAL  VALUES  ' 


''  Lambda  =  ",  E15.6E3)'3  RLAM 

"  b  =  " ,  ElS.Sea) ')  BVAL 

"  Beta  =  ",  E15.6E3J')  YROOr 


GOTO  999 
ERROR 

CALL  ERRMESSC lERR) 

CLOSE! 7) 

END 


{!■  {!■  -S' -s  ii- ->  *  *  -s-  if  •if  a  imifit »  a <■  -s  <■  if  if  i>  -s  if  if  if  if  if  <;•  a  -s-  if  •»  tf  -s  if  if  a  it  '>  *  <-  if  if  *  -s  if  if  -s- 

This  subroutine  reads  histogram  data  and  calculates  relative  freq 
SUBROUTINE  HISTORDC  FILENAME , DELTAA) 


implicit  double  precision  (A-H,0-Z) 

COMMON  /GGG/  NCL , FRELC 22S) , RHO( 22S) , XCl 2iS J , SUMA , SUMS , XBaR 
COMMON  /SSS/  OSHIFT,  NTOT,  ALPHA,  J20R3 

DIMENSION  IFABt  120) 

CHARACrER*30  FILENAME 

WHII E( * ,  'C A\)  ')  '  Enter  the  input  file  name 

READ!  ,  'C  A)  '  )  FILENAME 

OPEN! 6,  KILE  =  FILENAME,  STATUS  =  'DID') 

Read  absolute  frequencies 
READ! 8 ,  'C  A)  ')  COMMENT 
READC0,*)  A,  DELTAA 
READ! 8,  NUMCLASS 

BEACC  8,  =')  f  IFABC  ICNT)  ,  ICNT=1  .  MIJMC- AS.S'; 

NCL  =  NUMCLASS 

NTOT  =  0 

DO  100  JCNT  -  1,  NUMCLASS 

NTOT  -  NTOT  +  IFABC JCNT) 

CONTINUE 

DO  200  JCNT  *  1,  NUMCLASS 

FBELCJCNT)  =  DBLEl 1FA6C JCNT) }  /  uHLECNTOT) 

CONTINUE 


32 


a  n 


c 


300 


C 


200 


C 

C 

c 


Calculate  the  XC 

OSHIFT  -  ((-  DELTAA)/  2.0D0J  -  A 

DO  300  JCNT=  1  , NUMCLASS 

XC(JCNT)  -  DBLE(JCNT-I)  «•  DELTAA 

CONTINUE 

BETUHN 

END 


■-'•  ERRMEBS  writes  error  messeges  to  the  terminal 


SUBROUTINE  £BRMESS( lERR) 

IF  ( lERH . EQ . 0)  GOTO  500 
WRITE(*.<i-3  '  TERMINAL  ERROR!!!' 

IF  ( lEHH.EQ.  1j  THEN 

WRITE(*,")  '  No  7erocroS3ing  found  for  W0  ' 
EL3EIF  (IERfl.EQ.2j  THEN 

WRI TEC  '  Root  of  h  not  bracKetted' 

ELSEIF  (IERR.EQ.3)  THEN 

WRITEC '  h  does  not  converge  tc  a  zero' 
ELSEIF  (IEflR.EQ.4)  THEN 

WRITEC^f,*)  '  g-h  root  intersection  not  found' 

ENOIF 

RETURN 

END 


■•s  GFUNCf  X,  beta)  Calculates  the  derivative  g(c,beta) 

*  1  Sept  88 

*  mod  2  Dec  68 

Res.  Dir.,  C.  E.  Hall,  Jr. 

. .{(I  ^ j  jj.  i'f  ^  ^  ^  ^  ^  j-.  >«.  .;fc  <«.  .a  .y.  j.  tj. ^  ^  ^  ^  -a  ^  ^ 


i-'UNCTxUN  GFUNCC  X  ,  BETA) 

IMHLICIT  DOUBLE  PRECISION  CA-H,0-Z) 

COMMON  /GGG/  NCL , FHELi 225) , RHOC 225J  , XC( 228) , SUMA, SUMB , XBAR 
COMMON  /SS3/  DSHIFT,  NTOT,  ALPHA,  J20R3 

0  H  1 F  1  =  X 

CALL  SUMABC XSHIFT , BETA) 

3CETA  =  DLOG(BETA)  -  PSICBETA) 

GFUNC  -  (  OBLEI,  NTOT)  »X8AR-BUMB)  *BETA'■^SBETA  - 
(  uBLEC  NTOT)  <fXBAB  -  OBLEC  N  TOT  )  ••  XSHIF  T ) 


RETURN 

C.ND 
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CJ  CJ  o 


C  »«  e  -tm-  w-c;-  *  vHJ-o  •»**«•  if  ->»»»*  v-  ■!.-  a-  -nnnnnt  o  -e-'j-  *  ij-  <1  *  ■»  ::-  •>  it  «• »  *  a  *  * 

*  HFIJNC(BETA)  Calculates  the  derivative  hCc.beta) 

»  1  Sept  88 

*  13  Oct  88  mod  for  three  function  root  finding 

C  *  Res.  Dir.,  C.  E.  Hall,  Jr. 

C  S-it-tt-iUt-tHHf-a- a  if*  O*  it  if  iHr  if  '■•if  ii-fH!-*  it a  ****** 


FUNCTION  HFUNCC X, BETA) 

IMPLICIT  DOUBLE  PRECISION  CA-H,0-2) 

COMMON  /GGG/  NCL , FRELl 225)  , RHO(  225)  , XC(  225 J  , SUMA , SUMS ,  XOAL 
COMMON  /SSS/  DSHIFT,  NTOT,  ALPHA,  J20H3 

XSHIFT  =■  X 

IF  CJ20H3.EQ.3)  THEN 

SBETA  =  DLQGCBETA)  -  PSI(8ETA) 

YLOG  =  SBETA  /  (XSHIFT  -  XSAHJ 

SUM  =  a.oa 

00  100  JCNT-1 , NUL 

SUM  -  SUM  +  FREL(  JCNT)  xdEXPI  YLOG'-i  XCt  JCNT) -XC(  1  )  )  ) 
100  CONTINUE 

HFUNC  -  0EXP( YL0G»( XSHIFT-XC( 1 ) ) )  -  SUM 

ELSEIF  (J20R3.EQ.2)  THEN 

BVAL  -  BETA  *  ALPHA 
TEMPI  =  0.00 
TEMP2  '  0.00 

DO  200  JCNT  =  1 , NCL 

TERM  =  f-REL(JCNT)  /  DEXP(  XC(  JCNT )/ BVAL ) 

TEMPI  «  TEMPI  +  TERM 

T£UP2  =  TEMP2  +  TERM  XC(  JCNT) 

200  CONTINUE 

HFUNC  =  XBAR  -  BVAL  -  TEMP2  /  TEMPI 

ENOIF 

HE  TUHN 
END 


C  SUMA8(  X ,  BE  T  A )  Caicuiatas  the  series  sums  A  and  0 

C  1  S  e  p  t  0  8 

C  '  Res, Dir., C.E. Hall, Jr. 

C  *  it  if  if  if  if  a  *  *  if  if  if  if  if  if  if  if  if  if  -a  if  if  if  if  if  if  if  if  if  if  if  if  if  if  if  if  i.  if  if  if  if  *  if  if  '-f  if  *  *  'f  if  if  if  if  if  f  'f  if  *  if  if  if  if 


SUBROUTINE  SUMAB( X.BETA) 

implicit  double  precision  (A-H.O-Z) 
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COMMON  /GGG/  NCL , FREL( 225) , RHO( 225) , XCC 225 J , SUMA , SUMS , XBAR 
COMMON  /SSS/  DSHIFT,  NTQT,  ALPHA,  J20R3 

X6HIFT  -  X 

SBETA  -  DLOG(BETA)  -  PSICBETA) 

SUMA  =  0.00 
SUMS  -  0.00 

DO  410  J  -  1 , NCL 

SUMA  "  SUMA  4-  DBLEC  NTOT)  <-‘FREL(  J)  DEXP(  (  (XSHIFT  -  XCIJJ  ) 
CXBAR  -  XSHIFT)  )  '-f  SBETA  ) 

BUMS  -  SUMB  +  DaLE(NTOT)*EBEL( J)  *  XC( J)  -  DEXP( 

C  (XSHIFT  -  XC(J))  /  (X8AR  -  XSHIFT)  )  *  SBETA  ) 


410  CONTINUE 

RETURN 

ENO 

Q  ^  .{L  X}.  ^  ^  ^  ^  ^  ^  ^  M  ^  ^  .  V  {j.  {A  ^ 

PSI2  CALCULATES  VALUE  OF  PSI  FUNTION  FOR  GIVEN  ARGUMENT 

^  J'.  .K  ^  JL  .^4.  ^  ^  .{L  J}  .>4  .X.  ^  .>v  ^  ^  ^  M  .g.  ^  .n.  ^  ^  ^  ^  ^  M 

FUNCTION  PSI( YTX) 

IMPLICIT  DOUBLE  PHECISIQN( A-H , 0-Z) 

DATA  U,T1,T2,T3,T4,TS/1O0,1 2D0, 1000,2100, 2000 ,1100/ 

H=0 . 00 

XTX  =  YTX 

goto  200 

100  R=R+U/XTX 

XTX=XTX+U 

200  IF  (XTX. LE. 201)  goto  100 

Q  =  U/  (  XTXifXTX) 

PS1  =  LJ--H  Q'-H  Q^U  -Q/TS  +  U/T4)  -U/T3)  +U/T2) 

PSI=DLOG( XTX) -H-U/( EOO^XTX) +( PSI-U) *Q/T1 

RETURN 

ENO 

^  vj,  ^  xj.  .ij  jj.  ^  J/f  J}.  >{.  JL  >  •»>  i/e  'if  is-  i'r  '3-  *  if  ii  H  <K-  <-  if  ii  -ti-  -XJ  'Tr  if  i>  if  ->  -y*  ii-  it  it  H-  it 

SUBROUTINE  SETLIW(XMIN,XCEN,  XMAX  ,  DELxPOS  ,  OELXNEG.) 

IMPLICIT  DOUBLE  PRECISION  (A-H, 0-Z) 

COMMON  /SSS/  DSHIFT,  NTOT,  ALPHA,  J20R3 

COMMON  /GGG/  NCL  ,  FREL(  225)  , RH0(  225)  , XC(  225)  , SUMA , SUMB . XbAH 

XMIN  -  0.00 
XMAX  =  XBAR 

XCEN  =  XMAX  -  0.500  «  (  XC(2)-XCv1)  ) 
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nnno  ro  onnn  -*  n  oonono 


0ELXP08 

DELXNEG 


XMAX 

XCEN 


XCEN 

XMIN 


RETURN 

END 

*■  CALXBAR(  .]  calculatea  XBAR  for  the  gumbel 

*  distribution 

*  2  Dec  88 

*  C  E  Hall,  Jr,  Research  Directorate,  HOtC,  MICOM 

-;,m:  if  if  if  if  if if  *  'ri*  -*  .,■  ;f  ;.'  -;t  ,>  ;.- 


SUBROUTINE  CALXBAB 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-2) 

COMMON  /GGG/  NCL , EREL(  225)  , RHO(  22S)  , XCl  225)  , SUMA , SUMB , XBAR 
COMMON  /SSS/  DSHIFT,  NTOT,  ALPHA,  J20R3 

Calculate  the  moments 
XBAR  =0.00 
DO  100  ICNT=1 ,NCL 

XBAR  -  XBAR  +  PRELl ICNT)  «  XC( ICNT) 

00  CONTINUE 

RETURN 

END 

C  -If -If -tnnnnnnK^-if-lnnnfiHf -if  it  ■> -tut  innf  -If -itinnt-int -If -If  ■«■ 

<>■  R0QTG(XH0QT,YR00T,IERR) 

•»  31  Aug  88 

Res.  Dir.,  C  E  Hail,  Jr. 

^  >4.  ^  ^  .M.  ^  ^  .v.  ^  ^  .w.  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  {«,  «>y  ^  ^  ^  ^ 

SUBROUTINE  ROOTGC  XROOl  , YROOT , lERR) 

implicit  DOUBLE  PRECISION  (A-H,0-2) 

CALL  RIGBRACKC  XA, YA,XB, YB, I ERR, I0FLAG) 

IF  ( I0FLAG . NE . 0)  GOTO  100 
IF  ( IERn.NE.0  )  GOTO  200 

CALL  GIS4BTG(  XA, YA,xa, YB.XROOT, YROOT, lERR) 

GOTO  200 

00  ;<H(loi  =  XA 

YROOT  "  YA 

00  RETURN 

END 


•Or 


Jt  V J  ^  .Jj.  .JJ.  JJ,  .J J  ^  ^  ,K  JJ.  ^  .J}.  .JJ,  ^  ^  iJJ.  ^  {J.  ^  ^  ^  {J.  /b  ^  ^  ^  ,  L  ^  ^  X 

C VALUER  JCNT, 0ELXP0S,0ELXNEG, XPOS, XNEG) 

This  subroutine  calculates  the  positive  and  negative 
shift  values. 
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C  *  31  Aug  08 

*  Has.  Dir.,  C.  E.  Hall,  Jr. 

•It  *  -sum- »*■»*<!••»»«■» -IH!- «*•»****  iHJ- ■» 

SUBROUTINE  CVALUEC  JCNT  ,  OELXPOS  ,  OEL.XNEG ,  XPOB  ,  XNEG  ,  XCEN) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-2) 

TERM  -  DEXP(  DBLE(JCNT)  *OLOG(2.0D0)  ) 

XPOS  -  XCEN  +  TERM  '>  OELXPOS  /  (  TERM  +  99 . 0D0  ) 

XNEG  =  XCEN  -  TERM  DELXNEG  /  (  TERM  +  99 . 0D0  ) 

RETURN 
END 


C  *  RTGBRACK(.) 

C  *  This  subroutine  brackets  the  root  of  the  function  g 

C  #  31  Aug  88 

C  *  Res.  Dir.,  C.  E.  Hall,  Jr. 

C  ■it'C--it-it-in,‘it-itit->jti;--!.>-'Stit-!}'!t  It -a- -a- «•»*•>•«■«■»•»•**■»■•»«•«• -It -if  •? -it  It  it  J  ■it -It -it  it  it 


SUBROUTINE  RTGBRACK( XA , Y A , XB , YB , lERR , I0FLAG) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

JCNT  =  0 
I0FLAG  -  0 
IRLFLG  =  0 
MAXJ  "  30 

CALL  SETLIMC  XMIN , XCEN , XMAX , OELXPOS , DELXNEG) 

XC  =  XCEN 

CALL  ROOTHC XC, YC, lERR) 

IF( lERR . NE . 0)  GOTO  400 
ZC  =  GFUNC(XC,YC) 

IF  (ZC.EQ.0.00)  GOTO  390 

C  Calculation  loop 

200  JCNT  -  JCNT  -t  1 

CALL  CVALUEC  JCNT , OELXPOS , DELXNEG , XPOS , XNEG , XCEN) 

C  Check  positive  c  shift 

IF  ( C IHLFLG.EQ.2) .OR.C IRLFLG.EQ.3))  GOTO  235 
CALL  fiOOTHf XPOS , YP , lERRJ 
IF  (  lERR  .  NE . 0)  THEN 
lERR  =  0 

IRLFLG  -  IHLFLG  +  2 
GOTO  225 

ENDIF 

ZP  -  GFUNCC XPOS , YP) 

IF  (ZP.EQ.0.00)  GOTO  390 

IF  (  (  ZP^^ZC)  .  LT  .  0 . 00)  GOTO  350 
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c 

22S 


250 


C 

350 


C 

360 


C 

3'j0 


C 

395 


4  03 


C 

c 

c 

c 

c 

c 


Check  negative  c  shift 

IF  ( ( IRLFLG.EQ. 1) .0R.( IflLFLG.EQ.3) )  GOTO  250 
CALL  ROOTHC XNEG, YN, lERR) 

IF  ( lERR . NE . 0)  THEN 

TCRR  .  PI 

IRLFLG  =  IRLFLG  +  1 
GOTO  250 

ENDIF 

ZN  =  GFUNCC XNEG, YN) 

IF  (ZN.EQ.0.00)  GOTO  395 

IF  (  t  ZN-:ZC)  .  LT  .  0  .  D0)  GOTO  360 

IF  ((.  IRLFLG.  NE  .  3)  .ANO  .(  JCNT  .LT  .MAXJ)  )  GOTO  200 

lERR  =  4 
GOTO  400 


Sign  change  in  positive  J-1,J  bracket 
XB  =  XROS 
YE  =  YP 

J.WIN1  =  JCNT  -1 

CALL  CVALUEIJMIN 1 , DELXROS , DELXNEG , XA , XNEG , XCEN) 
CALL  ROOTHC XA, YA, lERR) 

GOTO  400 


Sign  change  in  negative  j,j-1  bracket 
XA  =  XNEG 
YA  =  YN 

JMIN1  =  JCNT  -  1 

CALL  CVALUEC  JMIN1 , OELXPOS , DELXNEG , XROS , XB , XCEN) 
CALL  ROOTHC X8, YB, lERR) 

GOTO  400 

Zero  found 
I0FLAG  »  1 
XA  =  XPOS 
fA  =  YP 
GOTO  400 

Zero  founo 
I0FLAG  =  1 
XA  =  XNEG 
YA  =  YN 

R  c.  T  C  R  T'J 
END 


■K-  if  ■!>  *  «•  *  if  i;-  -:f  -U  V.-  if  if  if  if  -f  if  if  if  if  if  if  if  if  if  if  'f  if  if  it  if  if  if  if  if  if  if  if  if  if  if  it  it  -f  i.-  ii-  if  if 

BIS4RTG  -  Bisects  the  interval  until  delta  X,  and 
delta  V  are  both  less  chan  1 .0E-7 
31  Aug  68 

Res.  Oir.,  C.  E.  Hall,  Jr. 


if  if  i.' 
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non 


SUBROUTINE  BIS4BTG(  XA , Y A , XB , YB , XM , YM , I ERR j 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

TOL  =  1 .00-07 

ZA  =  GPUNC(XA,YA) 

ZB  =  GFUNC[XB,YB) 

Bisection  loop 

XC  =  ( XA  +  XB)  /  2 . 000 

CALL  HOOTHC XC, YC, lERR) 

IF  ( IERR.NE.0)  GOTO  300 

ZC  =  GFUNC(XC,YC) 

IF  (ZC.EQ.0.D0)  GOTO  250 

IF  (  (  ZC<>ZA)  .  GT  .0 . 00)  THEN 
XA  =  XC 
YA  “  YC 
ZA  =  ZC 

ELSE 

XB  =■  XC 
YB  =  YC 
ZB  =  ZC 

ENOIF 

OELTAX  =  DABS(XB-XA) 

OELTAY  =  OABS(YB-YA) 

IF!  ( OELTAX. GE. TOL)  .OR.  ( OELTAY . G£ . TUL)  )  GOTO  100 

XM  =  [  XA  XB)  /  2.00 

YM  =  (.YA  +  Ya)  /  2.O0 

GOTO  300 

Exact  zero  found 
XM  ■=  XC 

YM  =  YC 

RETURN 

END 

^  .'/f  j/f  Jj.  .»<.  Jj.  J/.  zv  j>/.  .J,  z}.  ^  ZL  ZV  zv  ^  zj.  zj. .^j.  zj  z'. 

"■  RTHORACK  -  Root  H  Brocketing  oubrootiiTe 

This  subroutine  brackets  the  root  of  only  i;re  R  f 
3  1  A  u  g  8  8 

Res.  Dir.,  Dr.  C.  E.  Hall,  Jr. 

;;  •>  «■  i-  *  <•  *  *  *  *  c- » -3-  in;-  *  im-  o-  *  3  3333  3  c  *  3  -3  3  -3  3  -3  -3  3  3  3  *  3  3- 


SUBROUTINE  RTHBRACKC  XV,YA,ZA,YB,ZB,IERHJ 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 


Set  starting  parameters 


39 


n  n 


Y0 

-  0.100 

DUPFAC  -  1  . 

61600 

MAXITEH  -  1 

9 

lEHR  -  0 

Set 

initial 

values 

YA 

=  Y0 

ZA 

=  HFUNCC 

XV, YA) 

ITER  =  0 

Eva 

1 uat ion 

1  oop 

ITER  =  ITER 

+  1 

YB 

=  YA  *  DUPFAC 

ZB 

=•  HFUNCC 

XV, YB) 

IF 

(  ZB.EQ.0 

.00)  GOTO 

150 

IF 

(  (  ZA^»ZB) 

. LT . 0 . 00) 

GOTO 

YA 

»  YB 

ZA 

=  ZB 

IF 

C ITER . LT 

.MAXI  TER) 

GOTO 

MAX 

ITER  exc 

eeded  set 

error 

lERR  =  4 

ISa  RETURN 

END 


C 


i; 


•?:*  *  •»>  *  >«•  <5*  ♦>  -U-  '.J-  <»>  <«■  -S’  -.J 


'>-«?■?«  -» *•  <:•  <>•  •«•  •:>  •> 


BOOTH  -  Finds  the  root  of  h 
31  Aug  88 

Res.  Oir.,  Or.  C.  E,  Hail,  Jr. 

*  ^  a-  •;:•  ii  -j}-  •>  •{<•  •:>  <»>  <:•  <;•  -ji-  •>  <:•  -jj-  -u*  a  •>  ^r-  -jj  « *{:•  <:•  •::•  x- 


v*.  .K.  M  ^  ,•;.  oj.  ^  .'j  ,;;. 


SUBROUTINE  HOO TH(.  X V  ,  Y ROO  r  ,  lERR) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

CALL  HTHBRACKC  XV , Y A , ZA , YB , ZB , lERR; 

IF( lERR . NE . 0)  GOTO  200 

CALL  BROOTC  XV , YA , ZA , YB , ZB , YROOT , lERR) 

200  RETURN 

EMC 


SUBROUTINE  BROOTf  XV , Y A , ZA , YB , ZB , YROOT , lERR) 
implicit  DOUBLE  PRECISION  ( P-Z) 

DOUBLE  PRECISION  EPSI.HFUNC 
PARAMETER  (  M  AX  I  T=- 1  00  ,  EPS  1 -3 . 0-6) 

ZC  =  ZB 

00  610  NIT=1 .MAXIT 

IP  ( ZB*ZC . GT . 000)  THEN 
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YC  =  YA 
ZC-ZA 
YD=YB-YA 
YE-YD 
ENDIF 

IF(DABS(ZC) .LT.DA8SCZB))  THEN 
YA-YB 
YB  =  YC 
YC-YA 
ZA  =  ZB 
ZB  =  ZC 
ZC  =  ZA 
ENDIF 

TOL1-=2.O0*EPSl'»DAaSC  YB)  +  .SD0*EPSI 
YM=.5D0*(  YC-YB) 

IF  { DABS( YM) . LE . TOL 1 . OR . ZB . EQ . 0D0)  THEN 
YROOT-YB 
RETURN 
ENDIF 

IF( DAaS(  YE)  . GE . TOL1 . AND . DABS( ZA)  .GT . DABSf  ZB) ) 
s=za/ZA 

IFCYA.EQ.YC)  then 
T=*2D0<fYM*S 
Q= 1D0-S 
ELSE 

Q=ZA/ZC 

R=ZB/ZC 

T=S^^(  2O0«YM*Q'-»(  Q-R)  -(  YB-YA)  »(  R-1  .D0)  ) 

Q  =  (Q-1.D0)*(R-1.D0)«'(S-1. 00) 

ENDIF 

IF( T.GT.0.D0)  Q«-Q 
T=DABS( T) 

Y3=3O0*YM'-"Q-DA8S(  T0L1-"-Q) 

IF(  2  .  O0"‘T  .  LT  .  0MIN1(  Y3,  OABS(  YE»0)  )  )  THEN 
YE  =  YO 
YD=T/Q 
ELSE 
YD  =  YM 
YE=YD 
ENDIF 
ELSE 
YD  =  YM 
YE  =  YO 
ENDIF 
YA  =  YB 
ZA  =  ZH 

IF( DABSC  YD  I  . GT . TOL  1  )  THEN 
YB=Ya+YD 
ELSE 

YB-YB-t-DSIGN(  TOLl  ,  YM) 

ENDIF 

ZB-HFUNC( XV , YB) 

510  CONTINUE 
IEHR=S 
RETURN 
END 


THEN 
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c 

c 

c 

c 

c 

c 

c 

c 

c 


■»  *  -it  *  *  **■!!■•!!•  -a-  -u-  innt  -iHi-  -ii- » it  -K-  *  *  s  -K-  it  ■»  •!■  «•  •;:•  <;■  ■  *  if  if  if  if  if  -r-  if  if  if  if  if  if  if  if  if  it 

*■  This  program  calculates  the  gumbel  pdf  from  the 

gumbel  parameters 

*  powell's  method. 

*  Dr.  C.  E.  Hall,  Jr. 

*  Res.  Dir.,  ROEG,  MICOM 

*  17  Oct  1988 

*■  mod  5  Dec  1988 

*  if  *  *  if  if  it  ifif  ififif  ifif  ifiHf  *  *  if  if  if  if  if  if  *  if  if  if*  if  a- ifit  if  *  if  it  if  *  if  if- if  if  *  if  if  if  if  «  if  if  it  *  if  if  if  if  it  if  ii- 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
CHARACTER»30  FILENAME 

COMMON  /GGG/  NCL,  FRELC225),  RHOt 225; ,  XC( 225) 
COMMON  /SSS/  XSHIFT,  NTOT,  J20R3 

OPENf  7, FILE- 'PRN  ') 

WRITE!  '  3  Parameters' 

WRITE(-','{A'\)')  '  Enter  the  Lambda  -  ' 

HEAD!'--,*)  RLAM0 

WRITE!  '(  A\)  ')  '  Enter  b-  ' 

HEAD!  ■»■,*)  SIGMA0 

WRITE!  A\)  ')  '  Enter  Beta-' 

HEAD!*,*)  RMU0 


WHITE!  '  2  Parameters' 

WRITE! '^  '( A\)  ')  '  Enter  the  Lambda-' 

READ!  RLAM1 

WRITE! *,'! A\)  ')  '  Enter  b- ' 

READ! '-‘.-s:-)  SIGMA  1 

HMU1  =  1.00 

C  Set  parameter  for  unshift  XC's 

J20R3  =  2 
XSHIFT  =  0.00 

CALL  HISTOHD! FILENAME .OtLTAA) 

C  WRITE  HEADER  TO  PRINTER 

WRlTE!7,if)  '  GUMBEL  DioTHIBUTlON  lAlCULA  P  i  uNo  ' 

WRITE! 7,  '! A,  IX, A)  ')  '  Data  file  is  ’.FILENAME 

WRITE!  7,  '!  IX)  ') 

WHITE! 7,900) 

900  FORMAT  !1X,30X,'3  Para  E s t ima t e  '  , 7X ,  '  2  Para  Estimate') 

WRITE! 7,910) 

910  FORMAT  !  13X,  1Hv,5X,6Hf! abs)  ,5X,6Hf!cal)  ,bX,4HX»»2,6X, 

6Hf!  cal)  ,6X,4HX»if2) 

ISY1  =  0 
SY2  =  0.00 
SY3  =  0.00 
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DO  500  ICNT=1,NCL 


Y1  -  DBLE(NTOT)  *  FREL( ICNT) 

Y2  -  OBLEC  NTOT)  ii'PDFC  XC(  ICNT)  ,  RLAM0,SIGMA0,  HMU0)  *OELTAA 
Y3  =  DBLE(  NTOT) *POF(  XC(  ICNT)  ,RLAM1  , SIGMA  1  , HMU1 J ^OELYAA 
SY2  =■  SY2  +  Y2 
SY3  =  SY3  +  Y3 

DELI  =  CY1  -  Y2}  *  (Y1  -  Y2)  /  Y2 
DEL2  =  (  Y  1  -  Y3J  (  Y  1  -  Y3J  /  Y3 
lY 1  =  NINT( Y 1 J 
ISY1  =  1SY1  +  IY1 

WRITE(7,901)  ICNT,  I Y 1 , Y2 , DEL  1 , Y 3 , DEl2 

901  FORMAT  ( 12X, I3,7X, 13, 1X,4( 5X,F6 . 2) ) 

500  CONTINUE 

WRITE  (7,902]  ISY1,3Y2,SY3 

902  FORMA  T(  12X,  'Totals' 4X, 13, 6X,F6. 2,  16X,F6.2.1 
WRITEC  7,  '(  1H1]  '] 

CLOSEC  7] 

END 


Q  rfj.  ^  ^  ^  ^  .jj.  vj  ^  ^  .jt.  ^  ^  ^  .5.  ^  vj.  ^  ^ ^  ^  ^  ^  a  .y.  .fl.  ^  .3.  ^ ^  ^  ^  ^ 

C  This  subroutine  reads  histogram  data  and  calculates  relative  freq 

SUBROUTINE  HISTORO( F I LENAME , DELTAA) 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-2) 

COMMON  /GGG/  NCL  ,  F HEL(  225 )  ,  RHO(  225)  ,  XC(  22S.I 
COMMON  /SSS/  XSHIFT,  NTOT,  J20R3 
UlMENSlUN  IFAB( 120) 

CHARACTER^«-30  FILENAME 

WP1TE(-'','(A\)')  '  Enter  the  input  file  name  -' 

BE AOi;  A)  ')  FILENAME 

OPcN(b,  rIL.  E  =  HILENAME.  STAIuS  “  'OLD') 

C  Bead  absolute  frequencies 

BEADi.  B  ,  '(  A)  '3  COMMENT 
B£AD(8,'>3  A,  DELTAA 
BEADI  y, '-‘3  NUMCLASS 

READ'S,--)  (  I.-ABi;  ICNTJ  ,  I  CN  T 1  .  NUMCLASS) 

NCO  =-  NUMCLAES 

NTOT  -  0 

DO  100  JCNT  “  1,  NUMCLASS 

NTOT  -  NTOT  +  IFAB(JCNT) 

ICIl  CONTINUE 

00  200  JCNT  »  1 ,  NUMCLASS 

tHELlJCNTj  -  DULEl  lEABl  JCNT)  )  /  DBLE(NIOr) 

200  CONTINUE 
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n  o 


C  Calculate  the  XC 

IF  (J20H3.EQ.2)  THEN 

XSHIFT  =  (-1.D0)  »  XSHIFT 
00  250  ICNT  -  1 ,  NUMCLASS 

XC(  ICNT)  -A  +  (.  0aLE(  ICNT)  -0.500)  ^OELTAA  +  XSHIFT 

250  CONTINUE 

ELSE 

XSHIFT  -  ((  -  OELTAA)/  2.000)  -  A 

00  300  JCNT=1 .NUMCLASS 

XC(JCNT)  »  OaLE(JCNT-l)  OELTAA 

300  CONTINUE 

tNOIF 

RETURN 

END 


C  *  ii-  * a-  -a-  -s-  *  ---  *  i;- <■  *  «•  «•  M  '!•  ->  •> -y-  <!•  it  ->  *  tt  it  it  it  i.- it  it  it  it  a  it  it  it 

C  *  The  Gumbel  Probability  Oensity  function 

p  jj.  J1.  Jj.  .>1.  jj.  .y.  jj.  .y.  jJ.  Jj.  jj.  .5^  ^  .JL  .{j. .{}.  ^  ^  ^  -y.  ij;.  .h.  .jl  i\  ^  -  .{c  /J.  ^  j. 


FUNCTION  PDFf X.RLAMDA.aa.BETA) 

IMPLICIT  double  precision  tA-H,0-2) 

XI  =  (X  -  RLAMOA  )  /  BB 

FI  =  BETA  -SJ  DLQGCBETA)  -  DLOG(  BB)  -  GALOGl.  BETA) 
F2  ”  BETA  «•  (  OEXP(  -1.00  *  XI  )  +  XI  ) 

PDF  =  DEXP(  FI  -  F2  ) 

RETURN 

ENJ 


C 


•  1  0 


■>  it  it  it  it  ‘i  « it  it  it  *  «•  it  it  it  *  -t  it  it  a 

‘t  LN  GAMMA 

^  .^j.  jj.  ^  j;.  jj.  ^  ^  ^  ^  ^ 

FUMCTION  GALOGCQTEMP) 


•>  -x-  •;{•  if  •;>  -k-  *  •>  <-  it  «•  it  it  it  it  -it « it  -Jt  it  it  «•  it  it 

it  it  it  it  it  it  it  it  it  it  it  it  it  -.t  *•  it  it  it  it  it  it  it  it  it  <•  it  ‘It  i:  it  v  •>  it  it  it  it 


IMPLICIT  DOUBLE  PRECISION  (G,Q-Z) 
PARAMETER  (  U=  1  .  00  ,  XML  IM  =■  1  4 . 00) 
GALOG  =  999.9999 


IFi  QTEMP . EQ . 0 . 00)  RETURN 
X  --  UTEMP 
XM  Ju “U 

XMUL  =  XMUL  ••  X 
X  —  X  +  u 

IF(  X.LT.XMLIM)  GOTO  1  10 
0  -  X  *  X 

y=l  U/21  .  D0  +  (  C  -  1 .00)  /28  .  D0  +  B.O0/(  99.  DC^^O)  )  /Q)  /Q 
Y-l.5.D0  +  lU/(-6.D0)+YJ/a)/(  60  .O0»X)  -X*(  U-DLOG(  Xl  ) 
GALOG- Y*DLOG(  DSG3HT(  8  .O0«-OATAN(  U)  /X)  )  -OLOG(.  XMUL.) 
RETURN 
EriO 
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